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$> • Abstract 

We study the role of the unstable equilibrium points in the transfer of matter in 

^r) ■ a galaxy using the potential of a rotating triaxial system. In particular, we study 

r^ ', the neighbourhood of these points for energy levels and for main model parameters 

where the zero velocity curves just open and form a bottleneck in the region. For 

these energies, the transfer of matter from the inner to the outer parts and vice versa 

starts being possible. We study how the dynamics around the unstable equilibrium 

points is driven, by performing a partial normal form scheme and by computing the 

^ invariant manifolds of periodic orbits and quasi-periodic orbits using the reduced 

Hamiltonian. In particular, we compute some homoclinic and heteroclinic orbits 

playing a crucial role. Our results also show that in slow rotating and/or axisym- 

metric systems the hyperbolic character of the equilibrium points is cancelled, so 

that no transfer of matter is possible through the bottleneck. 
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1 Introduction 



This paper focuses on the study of the dynamics around the hyperbolic equi- 
librium points, L\ and L 2 , of a non-axisymmetric galactic potential, i.e. a po- 
tential whose principal axes on the (x, y) plane are different. In particular, we 
analyse the role of the invariant manifolds associated with the unstable peri- 
odic orbits and with invariant tori around L\ and L 2 in the large scale transfer 
of matter within the system. For this purpose, we use both semi-analytical and 
numerical techniques to compute the invariant manifolds and to study the role 
they play in the global morphology. This is a modern approach using dynam- 
ical systems which has successfully been applied in celestial mechanics and 
astrodynamics (e.g. |l|, |2|, y, U, |5(). A similar technique has already been used 
to explain the spiral arm and ring morphology in barred galaxies p, |7|, |8( . Pre- 
vious theories believe that spiral arms are density waves in a disc galaxy [9j. 
The density waves propagate from the centre towards the principal resonances 
of the galaxy, where they damp down [lOj ■ Other replenishment theories have 



been proposed, therefore, to obtain long-lived spirals (see [TJJ for a review). 
Using this innovative approach, we obtain outer rings and long-lived spirals. 
This paper is intended to explain the geometrical behaviour of such structures 
in a general manner. 

Such studies have many applications in galactic dynamics. Elliptical galaxies 
can be triaxial, (i.e. their principal axes can be all different from each other) 
while disc galaxies can contain several triaxial components, rotating or non- 
rotating, and with widely varying degree of non-axisymmetry, such as haloes, 
bulges, bars, or oval discs. In this paper, we consider a galactic model that 
describes a triaxial system. By studying its dynamics within a wide range 
of parameters, we will explain the dynamics of all the components and make 
links with their morphology. The logarithmic potential is a very suitable model 
for our purposes, because it has a simple expression, it is analytic, its series 
expansion can be calculated up to a high order, and it admits both semi- 
analytical and numerical treatments. 

We choose a reference frame such that the origin of coordinates coincides with 
the centre of the galaxy and the triaxial system is fixed, i.e. a reference frame 
rotating with the system. The idea of the present analysis is based on the fact 
that if we consider a range of energies for which the particles are confined in 
the inner region defined by the zero velocity curve, in our case, for energies 
lower or equal to that of the unstable equilibrium points (Ej(Li)), then there 
is no possible transfer of matter from the inner region of the galaxy to the 
outer region, or vice versa (see Fig. 1(a) and (b)). In this paper, we focus 
on energies slightly larger than that of the unstable equilibrium point, for 
which an opening in form of a bottleneck appears. We will hereafter refer to 
this aperture as bottleneck. Transfer of matter can be possible through this 
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Fig. 1. Location of the equilibrium points and the zero velocity curves, (a) Zero 
velocity curves at an energy level smaller than the energy of L\ and Li- Three 
regions are thus defined, namely an inner (solid dark grey), an outer (solid white) 
and the forbidden region (hatched light grey), and no transfer of matter between 
the inner and outer regions is possible, (b) Zero velocity curves at the energy level 
of the unstable equilibrium points L\ and L2. (c) Zero velocity curves at an energy 
level larger than the energy of L\ and Li- A bottleneck appears and matter can 
transit from the inner to the outer regions (solid white) and vice versa. 

bottleneck (see Fig. 1(c)), and we study the objects that drive the motion in 
this region. 

Therefore, our first goal is to study the neighbourhood of the unstable equi- 
librium points in the energy range where the bottleneck appears. This energy 
range is suitable to apply the normal form technique, which is known as the 
reduction to the centre manifold. As previously mentioned, it has been success- 
fully applied to celestial mechanics problems, in particular to the Restricted 
Three Body Problem, e.g. [ll,|3|,|2|]. The reduced Hamiltonian obtained from the 
reduction to the centre manifold process gives a good qualitative description 
of the phase space near the equilibrium points and uncouples the centre man- 
ifold from its hyperbolic behaviour. The procedure is similar to the Birkhoff 
normal form, usually used to study the stability properties around the central 
equilibrium point of the galaxy and to compute the families of periodic orbits 



around it, e.g. [12, ll3|, [14 



Our second goal is to study the behaviour of the hyperbolic invariant man- 
ifolds associated with the orbits contained in the centre manifold; these are 
essentially invariant manifolds of periodic orbits and invariant tori. As is well- 
known, stable and unstable invariant manifolds are dynamical features, that 
are responsible for the global dynamics in a dynamical system. In order to 
obtain the global picture of the transfer of mass, it is necessary to compute 
the invariant manifolds and to look for intersections between different parts, 
that is, obtaining the possible heteroclinic and homoclinic connections. These 
type of connections have been recently studied in the Restricted Three Body 
Problem and applied to obtaining transit and non-transit orbits in celestial 
mechanics |4|, |5| . 



In Section 2, we describe the characteristics of the galactic model. In Sec- 



tion 3, we present the equations of motion and we study the effect of the main 
parameters on the linear behaviour of the equilibrium points. In Section 4, we 
explain in detail the reduction to the centre manifold in the particular case 
of the logarithmic potential and perform a study of the practical convergence 
of the reduced Hamiltonian. In Section 5, we compute the invariant mani- 
folds associated with periodic orbits and quasi-periodic orbits and we study 
the role they play in the transfer of matter. We also perform an analysis on 
the variation of the free parameters of the potential. Finally, in Section 6 we 
conclude. 



2 Galactic model 



We consider a galactic model that fulfils the regularity requirements needed 
for the reduction to the centre manifold and that describes a triaxial system. 



We have selected the logarithmic potential [15j], which is analytical and has 
the expression 



1 / v 2 z 2 ' 

$(x,y,z) = -v 2 fog[R 2 + x 2 + ^ + — ). (1) 



The parameters p$ and g$ are non-dimensional constants that determine the 
shape of the potential, p$ being the planar and g$ the vertical axial ratios. 
Without loss of generality, we can consider the potential orientated such that 
< g$ < p$ < 1. On the (x,y) plane, higher values of p$ represent more 
circular systems, independent of q$. Analogously, on the (x,z) plane, higher 
values of g$ represent more circular systems, independent of p$. For a stan- 
dard triaxial figure, we choose p<$> = 0.75 and g$ = 0.65. The constants vq and 
Ro simply set the scale for the velocity and length, respectively. The velocity 
Vq corresponds to a circular asymptotic velocity at infinity and we consider 
Vq = 200 km s _1 , while when i? is different from zero, the singularity of the po- 
tential in the origin is avoided. We choose Rq = 14.14kpc. The pattern speed 



is set to Q = 5kms _1 kpc~ , characteristic for slow rotating systems [15 
The density associated with the potential is obtained via Poisson's equation 
V 2 $(x, y, z) = 4:iiGp(x, y, z). The terms in the Laplacian equation associated 
with the x and y components are always positive, while the term associated 
with the ^-component can become negative for large values of z. This intro- 
duces a restriction on the parameters p$ and g$ of the form g| > p%/(l +p%)- 
Throughout this paper we will use these standard values to fix the model and 
we will refer to it as "Model 1" . In the left panel of Fig 2, we show the isoden- 
sity contours of "Model 1" , and, in the right panel, the corresponding circular 
velocity curve, i.e. the velocity of a hypothetical star on a circular orbit. We 
note that we obtain elliptical isodensity contours and a flat rotation curve, 
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Fig. 2. Characteristics of the logarithmic potential. Left panel: Isodensity curves for 
"model 1" with parameters Q = 5kms _1 kpc _ ,i?o = 14.14 kpc, vq = 200 kms -1 , 
p$ = 0.75, and q$ = 0.65. Right panel: Corresponding rotation curve. 

these characteristics being appropriate for the study of triaxial systems. 



3 Equations of motion and study of the linear stability of the un- 
stable equilibrium points 



As is usual in galactic models, we consider a Hamiltonian formulation that, in 
our case, has a Hamiltonian function given by: 



H(x,y,z,p x ,p y ,p z 



^i.Pl+Pl+P 2 z ) + ^i. x ^y^ z ) -tt(xp y -yp x ) =Ej, 



where (x, y,z) are the coordinate positions, (p x ,p y ,p z ) are the associated con- 
jugate momenta in a rotating reference frame, and Q is the angular velocity of 
the galaxy in the inertial system (for a more detailed explanation see Binney 
& Tremaine 15|). The equations of motion in Hamiltonian coordinates are 



x = p x + Qy 


Px 


y = p y -ttx 


Py 


z = Pz 


Pz 



-® y - Qp x 

^*- y. ■> 



where $ x = ||, & y — |p, and § z = ff • The effective potential is defined as 
$ e g(x, y, z) = $(x,y,z) — Q(xp y — yp x ). The surface $ e flF = Ej is called the 
zero velocity surface, and its cut with the z = plane is the zero velocity 
curve. All regions in which $ e g > Ej are forbidden to a star, so we call them 
forbidden regions (see Fig. 1). 

As is well-known, the effective potential associated with the galactic model 
has five critical points, located in the xy-plane, where 
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dx 



dy 



dz 



0. 



Due to their similarity to the corresponding points in the Restricted Three 
Body Problem, they are often called Lagrangian points. L\ and L 2 are lo- 
cated on the x-axis symmetrically with respect to the centre, L 3 is placed at 
the origin of coordinates, and L 4 and L 5 are located on the y-axis symmetri- 
cally with respect to the centre. The linearised motion around L\ and L 2 is 
characterised by the superposition of a saddle behaviour in the xy-plane and 
two oscillations, one in-plane and one out-of-plane. Therefore, L\ and L 2 are 
linearly unstable points and their behaviour is known as saddle x centre x 
centre. They are usually called the hyperbolic points. The linearised motion 
around L 3 , L 4 , and L 5 is characterised by the superposition of three oscilla- 
tions, two in-plane and one out-of-plane. This behaviour is also known as a 
centre x centre x centre behaviour, so they are linearly stable and they are 
usually called the elliptic points. 



3. 1 Variation of the free parameters 



As previously mentioned, we concentrate on energy levels slightly larger than 
the energies of the equilibrium points L x and L 2 . At these energies there is a 
bottleneck in the region between the two zero velocity curves and we want to 
know the type of orbits that transit through it. 

The model parameters influence the global morphology, even for energy levels 
for which the bottleneck is open. Here we study the effect of the variation of 
the free parameters on the linear behaviour of the equilibrium points L\ and 
L 2 . For this purpose, we fix the values of Rq = 14.14kpc and Vq = 200 km s -1 , 
as in "model 1", and study the influence of Q, p$, and q$, separately. We 
make families of models in which only one of the parameters is varied, while 
the others are kept fixed. For each of the models in the family, we study 
the linear behaviour of the equilibrium points Li(L 2 ). We expand the effective 
potential, $ e g, around one of these points and we retain only first order terms. 
For Lx(L 2 ), the eigenvalues of the differential matrix corresponding to the 
planar motion are of the form: ±A, ±ui, and those corresponding to the 
vertical motion of the form ±z/i. Note that A, u and v are positive in general 
non-equal real numbers. Therefore, L\{L 2 ) is a linearly unstable saddle point, 
since the two real eigenvalues are related to the hyperbolic behaviour, while 
the purely imaginary are associated with the elliptic motion. 

We first study the effect of the pattern speed, Q. We vary the value of Q within 
the range f2 = 0.001 — lOkms^ 1 kpc -1 , the low values being characteristic for 
slow rotators and the large values, characteristic of fast rotators. In Fig. 3(a), 
we show how the modulus of the eigenvalues changes with the pattern speed. 
We observe that as the pattern speed decreases, the eigenvalues tend to zero, 
thus cancelling the hyperbolic character of Li(L 2 ) and removing the possibility 



of having transit orbits. The pattern speed is, on the other hand, related to the 
position of the equilibrium points through the expression Vt 2 = r^i d y J , 

where ri is the distance from the centre to L\ and $(r) is the potential on 
the equatorial plane. As Q decreases, r L increases, so the equilibrium points 
move farther out from the centre. In Romero-Gomez et al. [a, 13], we consider a 
barred galaxy model, whose hyperbolic equilibrium points are in the vicinity 
of the bar ends, and we relate the spiral arms and rings in barred galaxies to 
the invariant manifolds of the periodic orbits around the hyperbolic equilib- 
rium points. Elliptical galaxies, however, are considered to rotate very slowly 
or even not rotate 15|, ll6j. These triaxial systems have the equilibrium points 
too farther out and, according to our results, their hyperbolic character is 
cancelled. Thus, no transfer or escape of matter is possible, which is in agree- 
ment with observations because elliptical galaxies do not present spiral arms 
or rings. 



In Fig. 3(b) we show the effect of the planar shape parameter, p$, on the 
eigenvalues of the differential matrix. We consider models where we vary p$ 
in the range 0.6 — 0.99. We consider several models with different values of q$, 
limiting ourselves, for each value of p$, to models with a g$ value satisfying 
the restriction q\ > p|/(l +p%), in order to avoid negative densities in the 
z— axis mentioned before. The value of Q is kept fixed as in "Model 1". In 
all cases, as p$ increases, the eigenvalues tend to decrease. In particular, the 
real eigenvalue, A, tends to zero thus cancelling the hyperbolic character of the 
Li(L 2 ) points. However, note that, as expected, for a fixed value of g$ (vertical 
axial ratio), the corresponding vertical eigenvalue, u, does not change as we 
vary p$ (planar axial ratio). A similar behaviour can be observed in Fig. 3(c), 
where we study the effect of the vertical shape parameter, q$. The eigenvalues 
tend to decrease as the vertical axial ratio increases. However, for a fixed 
value of the planar axial ratio, p$, the values of A and uj remain constant. 
Therefore, we can conclude that only the p<$> shape parameter has an influence 
on the planar structure of the galaxy, while the parameter g$ will influence 
the vertical behaviour. 
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Fig. 3. Effect of the free parameters on the eigenvalues of the differential matrix: A 
(in red solid line and crosses), to (in green dashed line and stars) and 7 (in blue dot 
dashed line and open circles. The rest of the parameters are kept fixed as in "Model 
1". (a) Variation of the pattern speed, f2. (b) Variation of the planar axial ratio, 
p$. (c) Variation of the vertical axial ratio, g$. 



In conclusion, the hyperbolic character of L\ and L2 is, as expected, weaker 
for slowly rotating systems or for near-axisymmetric ones, e.g elliptical galax- 
ies. The hyperbolic character is strong and the bottleneck is open to global 
structures, in the opposite cases, e.g. bars. 



4 Reduction to the centre manifold 



The reduction to the centre manifold is a suitable tool for studying the neigh- 
bourhood of the hyperbolic equilibrium points. It was first introduced in 
Gomez et al. |l| for the Restricted Three Body Problem and further detailed 
in Gomez et al. [3j. The procedure can be summarised in two steps. First, 
the Hamiltonian function is expanded in power series around the equilibrium 
point. Then, a partial normal form scheme is applied in order to uncouple 
(up to a high order) the hyperbolic directions from the elliptic ones. Now the 
truncated Hamiltonian has an invariant manifold tangent to the elliptic di- 
rections of the linear part. The restriction to the invariant manifold tangent 
to the elliptic directions leads to a Hamiltonian system with two degrees of 
freedom and an elliptic equilibrium point at the origin. This restriction to the 
manifold is the so-called reduction to the centre manifold and the study of the 
dynamics of the reduced Hamiltonian gives a qualitative description of the 
phase space near the equilibrium point. 

We want to note here that we do not use the Birkhoff Normal Form (hereafter, 
BNF). We are interested in performing a local study of a hyperbolic equilib- 
rium point in the largest region possible. In that sense, the reduction to the 
centre manifold can provide better results in a larger neighbourhood than 
BNF. Furthermore, as we will see in Section 5, we are interested in comput- 
ing the invariant manifolds associated with periodic orbits and quasi-periodic 
orbits and the normal form described here provides as well initial conditions 
to compute both invariant objects. 

We proceed by first obtaining the quadratic real normal form of the second 
order term of the Hamiltonian function. This is achieved by performing two 
changes of coordinates. The first one consists of a translation plus an homoth- 
ecy in order to place the new origin of coordinates on the equilibrium point 
and to set the unit length equal to the distance from the equilibrium point 
to the origin. The second change of variables allows us to write the second 
order part of the Hamiltonian in its real normal form. We then perform a 
complexification of the variables, in order to obtain a Hamiltonian of diagonal 
form. Finally, the normal form of higher order terms is computed. In the fol- 
lowing subsections, we describe in detail the steps performed and in the last 
subsection we study the practical convergence of the reduced Hamiltonian. 



4-1 Translation and homothecy 

As previously mentioned, the first change of coordinates consists of a transla- 
tion plus an homothecy. The translation sets the new origin of coordinates on 
any of the symmetric hyperbolic equilibrium points and the homothecy sets 
the unit of distance as the distance from the origin of the old system to the po- 
sition of the equilibrium point. Furthermore, the change has to be symplectic, 
i.e. it must preserve the Hamiltonian equations in the new variables. 

If we write (x', y', z', p' x ,p' y , p' z ) as the Hamiltonian barycentric coordinates and 
(x,y,z,p x ,p y ,p z ) as the Hamiltonian ones around the hyperbolic equilibrium 
point, the symplectic change can be written as: 

x = (x 1 - xl.)/') y = <y - y Li )/l z = z'/i 

Px = lij>' x - Px L% ) Py = lij>'y- Py Li ) Pz = 1P' Z 



z' 



where Li = (xL i ,yL i ,0,p XL _,p yL _,0), i = 1,2 are the coordinates of the hyper- 



bolic equilibrium point L t and 7 = Jx L . + y Li - Expanding the logarithmic 
potential up to second order, the expression for the Hamiltonian in the new 
coordinates around Li, % — 1, 2, is 



H(x,y,z,p x ,p y ,p z ) = H + - — (p 2 x +p 2 y +p 2 z ) + -vfa 2 ( ( — - -j^j x 2 + 

l (j_ ML\ 2 kc Li y Li J_ 2 \ 

+ ^ Kp\ K*p% ) V K 2 p\ XV + Kq% " J ' 

-Q(xp y + yp x ) + ... 

where H = ±v%lnK + \{p 2 XL . + P 2 yL ) + ^(yL,Px Li -^uPy^ is a constant and 
K = Rl + x 2 L . + %. 



4-2 The quadratic real normal form 



The second change of coordinates will express the second order term of the 
Hamiltonian in its real normal form H 2 = Xxp x + \u>i(y 2 + Py) + \^2{z 2 +p 2 ). 
This is accomplished via a linear and symplectic change performed around 
the hyperbolic equilibrium point L\. The analysis for L 2 is analogous and the 
results are symmetric. 



For Li, we have 



and 7 = x^ 1 



\ ( ^ ) ~ i2 °' VLl= °' ' Pa;Ll = °' ' Ps/Ll = ° 




We note that the planar motion is uncoupled from the vertical motion, so 
we perform the reduction to the real normal form in two steps, that is, first 
regarding the in-plane motion and then the out-of-plane motion. In the 2D 
case, the differential matrix around L\ becomes 



M 






a b N 


-a 


6 


c 


0a 





— d —a i 



where 



a 



Y 



jL 



s- 2 *°, 



rr 1 



2„,2 



d 



n 2 7 



The matrix M has two real, ±A, and two purely imaginary, ±a>ii, eigenval- 
ues. The eigenvectors associated with the real eigenvalues have the following 
expression 



v\ 



2 t« 



2^«-f 



\2 „2 



2 „2 \ , A ( , A 2 a 2 

^a A + 5 l-c+ — - — 



(2) 



while the eigenvector associated with the imaginary eigenvalue is 



v\ = u + iv 



r — Z± — °L 
U b b 







+ i 



2f«i 







2t2-cji + it 



u b 



Considering V J V, where V = (v + \ u v_\ v), we obtain the suitable scaling 
factors d\ and d Wl that give the final symplectic change of variables regarding 
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the in-plane motion: 

v +x u V-X 



V 



r d\ 



Regarding the out-of-plane motion, the symplectic transformation is given by 



7 V^ 

Z Pz~> ~ Pz, 



u 2 7 

where uj\ = -^V > always. The final change of coordinates sets the Hamilto- 
nian around the equilibrium point L\ with coordinates that express the second 
order term in its real normal form. That is, after performing the second change 
of coordinates, the second order term of the Hamiltonian is: 

H2 = Xx Px + ^(y>+pl) + f (z'+pfj. (3) 



For subsequent computational purposes, it is desirable to work with complex 
variables. Then, a third change of variables is needed in order to write the 
second order term of the Hamiltonian in the complex diagonal normal form. 



4 ■ 3 Complexification 



We are interested in writing the second order Hamiltonian in a diagonal form, 
in order to simplify the homologic equation that we encounter in the nor- 
mal form process. We perform a complexification of the variables related 
to y and z. That is, we perform a symplectic change of coordinates from 
real coordinates (x,y,z,p x ,p y ,p z ) to complex coordinates of the normal form 
(QijQ2,Q3,Pi,P2,P3) complexifying the pairs (y,p y ) and (z,p z ) by means of: 

y-ipy z- ip z 

qi = x q 2 = 7=^ q 3 ~ 



Py ~W Pz~ iz 

Pi =Px P2= 1= P3 - 



V2 ^ V2 

v -iy Pz-i 

In this way, the second order Hamiltonian has the form: 

H 2 = Xqipi + iuiq 2 p2 + ^2<M>3, (4) 

where A, lo\ and uj 2 are real positive numbers. 

11 



4-4 Normal form of higher order terms 



We have now written the second order term of the Hamiltonian, H 2, in its real 
and in its complex normal form. The final step to obtain the centre manifold 
consists of removing some monomials in the expansion of the Hamiltonian so 
that the final Hamiltonian has an invariant manifold tangent to the elliptic 
directions of H2. We apply the method known as the Lie Series method (see 
[2I, |3j and references therein for a detailed description in the Restricted Three 
Body Problem). Until the end of this section, we use the following notation. 
If x = (xi, . . . , x n ) is a vector of complex numbers and k = (fc l5 . . . , k n ) is an 
integer vector, we denote x k the term x± . . . x^ 1 and in this context we define 
0° as 1. We define \k\ as Y^j \kj\- 

The initial Hamiltonian is expanded around the hyperbolic equilibrium point 
in the complex coordinates for which the second order term is in diagonal form 
(Eq. 4). Thus, the Hamiltonian_J has the form 



#(q, p) = F 2 (q, p) + E #n(q, P 



n>3 

where H2 is given as in Eq. 4 and H n is a homogeneous polynomial of degree 
n of the form ^^hijq^piqfp^qfp^, where hij denotes hi u i 2t i 3t j lt j 2t j 3 . 

The Poisson bracket of two functions, F(q, p) and G(q, p), of position and 
momenta is defined as 



M .t(^-^i. 



'dF_d^_d£dG\ 

dqidpi dpidqj 



If F and G are homogeneous polynomials of degree r and s respectively, then 
{F, G} is a homogeneous polynomial of degree r + s — 2. If a function G(q, p) 
is considered a Hamiltonian then, the function H defined by 

H = H + {H, G} + i«tf, G}, G} + ^{{{H, G}, G}, G} + ..., 

is the result of applying a time one map canonical transformation under the 
flow of G to the Hamiltonian H. The Hamiltonian G is usually called the 
generating function. If we choose a homogeneous polynomial of degree 3, G3, 
as a generating function, then the homogeneous polynomials of degree n, H n , 
such that H = J2 n >2 H n are given by 



Bold letters denote vector notation. 
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H 2 = H 2 , 

H 3 = H 3 + {H 2 , G 3 }, 

Hi = H4 + {H 3 , G 3 } + ^{{H 2 , G 3 }, G 3 }, 

In order to remove all the terms of order three in the new Hamiltonian, i.e. 
to have H 3 = 0, we must choose a generating function, G 3 , such that it 
solves the homological equation {H 2 , G 3 } = —H 3 . This procedure can be 
used recurrently trying to find homogeneous polynomials G n to remove non- 
resonant terms of the Hamiltonian. 

In our case, we are interested in removing the instability associated with the 
hyperbolic character of the Hamiltonian H. We note that the second order 
term of the Hamiltonian provides the linear part of the Hamiltonian equations. 
Therefore, the instability is associated with the term \q\p\. For the linear 
approximation of the Hamiltonian equations, the centre part can be obtained 
by setting q\ — pi — 0. If we want the trajectory to remain tangent to this 
space (i.e. to have qi(t) = Pi(t) = for all t > 0), then we need <7i(0) = 
p\ (0) = when adding the nonlinear terms. Due to the autonomous character 
of the Hamiltonian system, we will obtain qi(t) = pi{t) = for all t > 0. 
We remember that the Hamiltonian equations associated with a Hamiltonian 
H(q, p) are 

OH . _dH 

Qi t; i Pi 7; • 

opi oqt 

n^3 

Therefore, we can obtain the required condition, gi(0) = pi(0) = when 
9i(0) = Pi(0) = 0, if in the expansion of the Hamiltonian H there are no 
monomials with exponents (1, 0,i 2 , j 2 , i 3 , j 3 ) and (0, 1, i 2 , j 2 , i 3 , j 3 ). Different 
changes of variables can be used for such purpose, the most common being 
the ones cancelling the terms with i\ 7^ ji [3] or cancelling only terms with 
i\ + ji = 1. This second choice cancels the minimum number of terms in the 
Hamiltonian and is the one we have chosen in this paper. 

Thus, after the final change of variables the Hamiltonian can be written in 
the form: 

_ #(q,P) = H N (q,p) + J R 7V (q,p), 

where Hn(q, p) is a polynomial of degree iV without terms of i% + ji = 1 and 
Rn is a remainder of order N + 1, which is very small near Li(L 2 ). 



In particular, 
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Finally, using the inverse change of variables of the complexification, the trun- 
cated Hamiltonian, H n can be expanded in real form and we obtain 



iY 



#Jv(q, P) = H 2 (q, p) + J2 H n{<\, P) 

n=3 

where the second order term H 2 (q, p) is as in Eq. 3. 
4-5 Practical convergence 



In order to check the practical convergence of the truncated series, we inte- 
grate the same initial condition on the centre manifold using both the reduced 
Hamiltonian up to order 15 and the Hamiltonian in barycentric coordinates. 
At each time step, we compare the two position vectors by computing the 
norm of the difference. When this becomes greater than a given tolerance, we 
plot the initial condition and the time step it has reached. In Fig. 4, we show 
the practical convergence plots for four different energy levels and a tolerance 
of e = 10 -6 for "model 1". In Fig. 5, we plot the energy levels for the same 
model but with a tolerance of e = 10~ 9 . The energy of the equilibrium point is 
Ej(Li) = 130055.178. In each panel, we consider an energy level higher than 
the previous one and initial conditions on the centre manifold. We note that, 
due to the energy confinement, the initial conditions valid for a given energy 
level lie inside the planar Lyapunov orbit. As the energy increases, the time, 
t, at which the differences become greater than the tolerance, decreases. This 
is reflected in Figs. 4 and 5 by the colour palette. In all panels, the colour 
palette ranges from t — to t — 6.5, the green colour being associated with 
large times and the red colour, to small times. In Fig. 4, where we use e = 10~ 6 
as tolerance, we note that in the top left panel, the dominating colour is green, 
while as we increase the energy, the colours tend to yellow and red. We also 
observe that, as expected, as the tolerance decreases, the time also decreases, 
so in the bottom right panel of Fig. 5 the region is practically all orange and 
red. 

Note that for a tolerance of 10~ 6 , the practical convergence of the normal 
form is good for a wide range of energies, as all panels in Fig. 4 are basically 
green. Even for a more restrictive tolerance, 10 -9 , the normal form gives good 
results for a considerably range of energies, see top panels of Fig. 5. Also note 
that the numerical test we use here integrates initial conditions using both the 
reduced and the initial Hamiltonian. Since the initial Hamiltonian contains the 
instability, that grows exponentially in time, the time estimations we show 
here underestimate the real values of t for which the reduced Hamiltonian 
gives accurate results. 
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Fig. 4. Study on the practical convergence of the reduction to the centre manifold 
for a tolerance of e = 10 -6 in "model 1". In all panels, we associate to each initial 
condition and the final integration time to a colour in the palette. The energy levels 
used are 130060.178 (top left), 130300.178 (top right), 130600.178 (bottom left) and 
130900.178 (bottom right). 
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Fig. 5. As in Fig. 4 for a tolerance of e = 10 9 
5 Invariant objects around L\jLi. Transfer of matter 



In the following, we apply the methodology introduced in the previous section 
to compute the invariant manifolds associated with both periodic orbits and 
invariant tori. We fix the values of the potential and study the neighbourhood 
of the hyperbolic equilibrium point L\. We first compute the invariant curves 
around the equilibrium point on the section given by X = {z = 0} and within 
a range of energies. We then compute the invariant manifolds of the planar 
and vertical Lyapunov orbits and of the invariant tori. We also study the 
role invariant manifolds of periodic orbits have in the transfer of matter in 
the galaxy by computing the homoclinic and heteroclinic connections between 
the periodic orbits of a given energy level. Finally, we study the variation of 
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the free parameters and check their influence on the presence of homoclinic 
and heteroclinic orbits, thus, on the global shape produced by the invariant 
manifolds. 



5.1 Invariant manifolds of periodic orbits and quasi-periodic orbits 



Throughout this section we fix the model parameters to those of the standard 
"model 1" and we compute the invariant objects using the Hamiltonian in 
the normal form up to order 15 in the reduced coordinates. For clarity of the 
representation, however, at each step we perform the backwards change of 
variables to plot the results in the initial frame using barycentric coordinates. 

We compute the invariant curves on the section X and the planar Lyapunov 
orbit surrounding the equilibrium point, see top panel of Fig. 6. To obtain the 
invariant curves, we take initial conditions on the galactic plane of the initial 
frame of reference and we plot the cuts with the section. Each invariant curve 
on the section represents an invariant torus, as is shown in the top and bot- 
tom left panels of Fig. 6. In the bottom left panel, we plot the invariant torus 
corresponding to the invariant curve on the section X (red solid line) in the 
original coordinates, although the numerical integration is performed using 
the reduced Hamiltonian. Since the reduced Hamiltonian has the hyperbolic 
component decoupled from the elliptic one, we can perform long time inte- 
grations without the instability. For convenience, however, we use the normal 
form change of variables at each integration step to display the orbits in the 
original coordinates. In the bottom right panel, we plot the vertical Lyapunov 
orbit. We observe that for the range of energies for which the normal form of 
the Hamiltonian is accurate, the section X of the invariant curves has always 
the same aspect, i.e., the only macroscopic objects we see are the Lyapunov 
orbits and invariant tori. 

We can also compute the hyperbolic invariant manifolds associated with the 
planar and vertical Lyapunov orbits and with the invariant tori. We first com- 
pute the stable and unstable invariant manifolds associated with the planar 
Lyapunov periodic orbit around the equilibrium point L 1; W* and W™ re- 
spectively, using the truncated Hamiltonian up to order 15. According to the 
normal form scheme used, the initial conditions approximating the unstable in- 
variant manifold of a planar Lyapunov orbit have the form (qi,pi, q2,P2, <?3, Pz) = 
(±e, 0, q2,P2, 0, 0) and the initial conditions for the stable invariant manifold 
are (qi,Pi, <?2,P2, Q3,P3) = (0, ±e, <?2,P2, 0, 0). Note that invariant manifolds are 
sets of asymptotic orbits that tend to and depart from the periodic orbit. In the 
left panel of Fig. 7, we plot W^ and W^ for the energy level Ej = 130155.178 
(Ej(Li) = 130055.178) and e = 10~ 5 . We also plot the zero velocity curve 
defining the forbidden region. Here we observe the characteristic saddle be- 
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Fig. 6. Top panel: Invariant curves (black dots) and planar Lyapunov orbit (green 
solid line) around the equilibrium point L\ for the logarithmic potential with pa- 
rameters as in "model 1", and energy level Ej = 130100.178. The red solid line 
corresponds to the invariant torus plotted in the bottom left panel. Bottom left 
panel: Invariant torus (red solid line) corresponding to the (red solid line) invariant 
curve on the section X. Bottom right panel: Vertical Lyapunov orbit (blue dotted 
line) of the same energy. 

haviour of the invariant manifolds due to the hyperbolic behaviour of the 
periodic orbit. Also note that the invariant manifolds connect the inner region 
with the outer region, in the sense that the existence of invariant manifolds 
of periodic orbits implies also the existence of transit orbits between the two 
regions delimited by the zero velocity curves. In the right panel of Fig. 7, we 
plot the unstable invariant manifold of an invariant torus for the energy level 
Ej = 130155.178. The initial conditions to compute the unstable and stable 
invariant manifolds of invariant tori are easily provided by the normal form we 
use. That is, (qi,pi,q2,P2,q3,P3) = (± e , 0, g 2 ,P2, ^3,^3) for the unstable mani- 
fold, (0, ±e, q2,P2, q3,P3) for the stable manifold and e = 1CT 5 . We superimpose 
it to the unstable invariant manifold of the planar Lyapunov orbit of the same 
energy and we observe that they describe approximately the same loci in the 
(x, y)-projection, the invariant manifold of the invariant torus being located 
inside the invariant manifold of the planar periodic orbit. 

In Fig. 8, we plot the unstable invariant manifolds of both the planar, 71, and 
vertical, Si, Lyapunov orbits of the same energy level (Ej = 130155.178). The 
initial conditions used to compute the unstable invariant manifold of the verti- 
cal Lyapunov orbit have the form (q±,pi, <?2, £>2, <?3, Pz) = (± e ? 0, 0, 0, q 3 , p 3 ) and 
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Fig. 7. Invariant manifolds of a periodic orbit (left panel) and of a quasi-periodic 
orbit (right panel) of the logarithmic model for the energy level of Ej = 130155.178. 
Left panel: In the centre of the plot, a black solid line shows the planar Lyapunov 
orbit around L\. We plot the two branches of the unstable invariant manifold (red 
dotted lines), the two branches of the stable invariant manifold (green dotted lines), 
and the the zero velocity curves (grey solid lines) defining the forbidden region. 
Right panel: Invariant manifolds of an invariant torus, Wj, (in green), and of a 
periodic orbit, W^ x (in red). In black solid line, the planar Lyapunov orbit around 
L\ and, in grey, the zero velocity curves. 

the initial conditions for the stable invariant manifold are {qi,Pi, Q.21P21 Q3, P3) 
= (0, ±e, 0, 0, g 3 ,p 3 ), with e = 10~ 5 . In the left panel, we plot the invariant 
manifolds on the projected galactic plane. We observe that they both describe 
essentially the same loci, the invariant manifold of the vertical Lyapunov orbit 
being located inside the invariant manifold of the planar Lyapunov orbit in 
the configuration space. In the right panel, we show the three dimensional 
view of the same manifolds. We observe that the invariant manifolds of the 
planar Lyapunov orbits drive the motion dominating over the invariant mani- 
folds of the vertical Lyapunov orbits. This is due to the fact that the unstable 
component lies within the plane, since the two real eigenvalues are associated 
with the in-plane motion. We can, therefore, conclude that the motion in the 
neighbourhood of L\ and L2 is basically dominated by the invariant manifolds 
of planar periodic orbits. 



5.2 The role of invariant manifolds in the transfer of matter 



As we have seen, the motion in a neighbourhood of a hyperbolic point is 
mainly driven by the invariant manifolds of planar periodic orbits, so in this 
section we restrict ourselves to the z = plane. The next step is to study 
the global behaviour of the invariant manifolds associated with planar pe- 
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Fig. 8. Two different views of the invariant manifolds of the planar Lyapunov peri- 
odic orbit, W" (in red), and of the vertical Lyapunov periodic orbit, W$ (in blue), 
of the logarithmic model for the energy level of Ej = 130155.178. Left panel: (x,y) 
plane. Right panel: 3D view (not to scale). The black solid line and the yellow dot- 
dashed line show the planar and the vertical Lyapunov orbits around L\ , 71 and b\ , 
respectively. 

riodic orbits and determine their role in global structures. In particular, we 
are interested in knowing whether the invariant manifolds can drive particles 
from the neighbourhood of L\ to the neighbourhood of L2 or even connect 
the neighbourhood of L\ with itself using long-time trajectories. For clarity, 
as we did in previous sections, we plot the Poincare maps and the orbits in 
the initial frame using barycentric coordinates. 



A particle will be transferred from the vicinity of one Lagrangian point to the 
vicinity of the symmetric point if it is "trapped" first by W" and then by W* . 
To study these type of transitions, we use Poincare surfaces of section, that 
is, we draw the crossings of trajectories through a particular plane or surface 
in phase space. Depending on the purposes of our study, some surfaces will 
be more suitable than others, but the methodology is the same. Let us take 
as an example the surface of section S defined by y = with x > 0; that is, 
we consider the orbits when they cut the plane y = having a positive value 
for the x coordinate. Let us consider this surface of section S for the stable 
and unstable invariant manifolds of a Lyapunov orbit around L2 (located in 
the x < side). Taking initial conditions on the manifold and integrating W" 
forward in time (resp. W? backwards in time) until the first encounter with 
S, we obtain the simple closed curves W^ 1 (resp. W^ 2 l ) that can be seen in 
Fig. 9a. Although the simple closed curves W^ 1 and W^ 1 are obtained as the 
natural result of intersecting the manifold tubes with a plane, it should be 
mentioned that further crossings (i.e. W™ 2 ,k and W^ 2 with k > 1) may well 
not have this simple structure, as can be seen in the Restricted Three Body 
Problem example provided by Gidea and Masdemont 17 . 



In the selected example, W^ 1 and W^ 1 are represented in (x,x) 



coordinates. 
It is important to note that a pair (x, x) in S defines an orbit in a unique way, 
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since y = and y is obtained from the energy level under study and in the sense 
of crossing S. By definition of invariant manifold, a point in W^nW^ 1 (black 
dots in Fig. 9a) represents a trajectory asymptotic to the Lyapunov orbit 7 
around L 2 both forward and backward in time and this trajectory is called 
a homoclinic orbit. In general, homoclinic orbits correspond to asymptotic 
trajectories, ip, such that ip e W" fl W* v i — 1,2. Thus, a homoclinic orbit 
departs asymptotically from the unstable Lyapunov periodic orbit 7 around 
Li and returns asymptotically to it, as can be seen in Fig. 9b. 
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a model with Q = 5 km s x kpc 
0.7 and g$ = 0.65 (Not to scale), (a) 



Fig. 9. Homoclinic 

v = 200kms~ 1 , R 

We plot the closed curves W 72 ' (red solid line) and W^ 2 (green dashed line) on the 

surface of section S defined by y = with x > 0. The black dots correspond to 

homoclinic orbits, (b) We plot the same invariant manifolds as in panel (a), now 

in the configuration space (x,y). The black curves correspond to the homoclinic 

orbits. 



Heteroclinic orbits, on the other hand, are defined as asymptotic trajectories, 
ip', such that ip' G W™. fl W*., i 7^ j, i,j = 1,2. Thus, a heteroclinic orbit 
departs asymptotically from the periodic orbit 7 around L, and asymptotically 
approaches the corresponding Lyapunov periodic orbit of the same energy 
around the Lagrangian point at the opposite end of the bar region Lj, i 7^ j; 
a suitable surface of section, S', for this computation can be the plane x = 
with y > 0. We hereafter refer to interior region as the ellipticahlike shape 
defined by the iso-effective potential curves encircling the centre and roughly 
passing through L\ and L 2 . In Fig. 10a, we plot the closed curves W^ 1 and 
W^ 1 on the surface S'. Note that there are two intersection points, therefore 
this model presents two heteroclinic orbits that will connect asymptotically 
one side of the interior region with the opposite, as can be seen in Fig. 10b. 



The existence of these type of connections determines the way matter is trans- 
ferred within the galaxy. By definition, if there are heteroclinic connections, 
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10. Heteroclinic 
200 km s" 1 , R = 



connections for a model with Q = 5 km s 1 kpc , 
14.14 kpc, pq, = 0.95 and q® = 0.85 (Not to scale), (a) 
We plot the closed curves W"/ (red solid line) and W^ (green dashed line) on the 
surface of section S' defined by x = with y > 0. The black dots correspond to 
heteroclinic orbits, (b) We plot the same invariant manifolds as in panel (a), now 
in the configuration space (x,y). The black curves correspond to the heteroclinic 
orbits. 

they create a tube, where particles can be trapped and travel from one side 
of the interior region to the opposite. The intersection of the tube with the 
section S' can be seen in Fig. 10a as the hatched area obtained from inter- 
secting the curves W^ 1 and PF? ' 1 . The global morphology in such models is 
reminiscent of that of Ri rings |l8l . |g, |3] . R\ rings are outer rings in barred 
galaxies whose semi-major axis is perpendicular to the one of the bar. They 
have a characteristic 9 or 8-shape (see Fig. 11a). If the model presents homo- 
clinic connections, particles can be trapped in tubes travelling from one end of 
the interior region to itself. The global morphology in this case is reminiscent 
of that of -R1-R2 rings 18|, 13] (see Fig. lie). R1R2 ringed galaxies present two 
outer rings, namely the Ri ring previously mentioned and an R 2 ring, whose 
semi-major axis is parallel to the bar semi-major axis (see Fig. lib). Note 
that if neither heteroclinic nor homoclinic connections exist, particles are not 
able to return to the vicinity of the interior region in an early stage. Thus the 
particles follow an escaping trajectory forming spiral arms [7J (see Fig. lid). 



In Sect. 5.1 we show that motion is mainly driven by the invariant manifolds of 
the planar Lyapunov orbits. Here we follow the study in Sect. 3.1 and we com- 
pute the invariant manifolds of families of models where we vary the pattern 
speed, Q, and the planar axial ratio, p$, i.e. the two dynamical parameters 
that have an effect on the potential in the galactic plane. The rest of the pa- 
rameters are fixed as in "model 1" . For each model, we compute the family of 
planar Lyapunov orbits in a range of energies for which they are unstable and 
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Fig. 11. (a), (b) and (c) schematic prototypes of outer rings from Buta & Crocker 
[18| reproduced by permission of the AAS. (d) schematic prototype of a spiral arm. 

the bottleneck through the zero velocity curve is opened. We compute the in- 
variant manifolds associated with the planar periodic orbits and we study the 
existence of possible heteroclinic and homoclinic connections. In each panel 
of Fig. 12, we plot the corresponding global morphology. In the left column 
we increase from top to bottom the value of the pattern speed. In all cases 
the global shape of the invariant manifolds is two spiral arms and the degree 
of openness is essentially independent of the value of the pattern speed. This 
can be seen in Fig. 13, where we plot a measure of the openness of the spiral, 
namely the ratio, R s , between the x-coordinate of the outer branch of the un- 
stable invariant manifold at the first cut with the y = axis and the absolute 
value of x-coordinate of the equilibrium point. In the right column of Fig. 12, 
we plot the invariant manifolds for models with different values of p$. When 
the system is strongly non-axisymmetric, no homoclinic or heteroclinic con- 
nections are present and the global structure is that of two spiral arms. When 
p$ = 0.7, the model presents homoclinic connections with both ends of the 
interior region. As previously mentioned, the global morphology reminds that 
of R1R2 rings. As p$ approaches to 1, that is, the system becomes prolate, 
the invariant manifolds tend to close, heteroclinic connections are present and 
the structure forms an R± ring. We stress the similarity of these patterns with 



the schematic prototypes of outer rings described by Buta & Crocker [18j. In 
Fig. 13, we also plot the degree of openness as the shape parameter p$ varies. 
As is illustrated in the right column of Fig. 12, the value of the ratio decreases 
as the value of p$ approaches to 1. 



6 Conclusions 



In this paper we use a suitable galactic potential (the logarithmic potential) 
to perform a semi-analytic study of the neighbourhood of the unstable equi- 
librium points L\ and L 2 . For a detailed study of their neighbourhood, we use 
a partial normal form scheme and we find that the main objects are the planar 
and vertical families of Lyapunov orbits and invariant tori. We compute the 
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Fig. 12. Effect of the pattern speed 0, (left panels) and of the planar shape parameter 
p$ (right panels) on the invariant manifolds. 

invariant manifolds associated with both periodic orbits and quasi-periodic or- 
bits using the reduced Hamiltonian and we find that the motion around L\ and 
L2 is mainly driven by the invariant manifolds of the planar Lyapunov orbits. 
However, we are also interested in the global structure of the galaxy. Thus, we 
study the possible homoclinic and heteroclinic connections between the pla- 
nar periodic orbits. For such a purpose, we use suitable Poincare surfaces of 
section. We note that this approach has successfully been used in celestial me- 
chanics and in this paper we apply it to a galactic dynamics problem, namely 
the formation of spiral arms and rings in barred galaxies, in comparison to 
other theories given so far. 

Here we are interested in the hyperbolic behaviour of L\ and L 2 and, partic- 
ularly, in determining the role the invariant manifolds play in the transfer of 
matter for energy levels where the zero velocity curves are open (i.e. a range 
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Fig. 13. Effect of the pattern speed (black solid line and labels in the bottom) and 
of the shape parameter p$ (red dashed line and labels in the top) on the degree of 
openness of the invariant manifolds, R s . 



of energies somewhat larger than the energy of L\ and L^) and a bottleneck 
appears around L\ and L 2 . For such purpose, we also study the influence of the 
main model parameters, namely the pattern speed, Q, and the shape param- 
eters, p$ and q$. The logarithmic models are suitable for describing triaxial 
systems such as haloes, bars, bulges in disc galaxies or elliptical galaxies. The 
structures we have constructed using the invariant manifolds, however, are 
not globally dependent on the model characteristics [fj, |7|. This implies, for 
example, that if a system had some degree of rotation, these kind of structures 
should be present. Elliptical galaxies are triaxial systems that do not present 
any external feature, i.e. they are an ellipsoidal distribution of matter with 
different degrees of ellipticity and they present neither spiral arms nor rings. 
Observations show that elliptical galaxies barely rotate or do not rotate at all 



as a figure [15l. Ilq|. Our results are in agreement with this statement. We have 
shown that models that rotate slowly or do not rotate cancel the hyperbolic 
behaviour of the equilibrium points and, thus, no transfer or escape of matter 
is possible. On the other hand, bars in disc galaxies are non-axisymmetric 
components usually characterised in the literature by elliptical distributions 



of density [19( although both observations and simulations show they might 



have more rectangular endings [20, 21]. It is well-known that bars rotate and 



observations show that spiral arms or rings emanate from the ends of the bar. 
This characteristic is also consistent with our results which show that if the 
system rotates at a given angular velocity, the hyperbolic equilibrium points 
are present, so that the invariant manifolds drive the motion, and therefore, 
set the global morphology to the galaxy. We have seen that, depending on the 
rotation velocity and the shape of the bar, the morphology will be that of a 
barred spiral galaxy or that of a barred ringed galaxy. 
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